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1. Introduction 

Recent development of algorithms and computational resources has enabled extensive studies 
of dynamical lattice QCD simulations at small quark masses where the chiral symmetry plays 
an essential role. Such studies are indispensable to investigate the chiral dynamics of QCD and 
to determine hadronic matrix elements with the precision required by present and future flavor 
physics. Among lattice fermion actions, the overlap fermion [jj], |||] has attractive features for these 
purposes. The overlap fermion operator is represented as 

D = m [l + Y 5 sign(H w (-m ))}, (1.1) 

where Hw = JsD w is the hermitian Wilson-Dirac operator with a large negative mass — iuq. This 
operator satisfies the Ginsparg- Wilson relation and thus holds an exact chiral symmetry on the 
lattice [Q, a\. The exact chiral symmetry significantly simplifies the structure of operator mixing in 
calculations of the matrix elements. The overlap fermion corresponds to the infinite N s limit of the 
domain-wall fermion, which means that one does not have to take care of the residual mass. On 
the other hand, numerical implementation of the overlap fermion is expensive, because of the eval- 
uation of the sign function of Hw- Furthermore, the discontinuity of the operator at A = 0, where 
X is an eigenvalue of Hw, makes the molecular-dynamics evolution much involved. Therefore, 
dynamical simulations of the overlap fermion have become feasible only recently with improved 
algorithms and large computational power. 

We are running a large-scale lattice QCD simulation project with 2 and 2+1 flavors of dynam- 
ical overlap fermions. Our physics goals are to explore the chiral regime of QCD with an exact 
chiral symmetry, and to compute hadronic matrix elements with controlled chiral extrapolation. 
Present simulations are performed with a spatial lattice size of 16 and a ~ 0.12 fm. For the gauge 
action, we adopt the Iwasaki's renormalization group improved action, as well as a topology fixing 
term which suppresses near-zero modes of Hw [§, 0, §1. Thus our simulation is performed in a 
fixed topological charge sector. By avoiding the discontinuity of the overlap operator at A = 0, the 
numerical cost of HMC is significantly reduced. 

This report explains our strategy of simulations and presents several recent results. In the next 
section, we address why the fixed topology simulation is feasible, and how it can extract physical 
observables with controlled systematic errors. Section || describes the numerical algorithms which 
are essential for dynamical overlap simulations. In Section |], our simulation set-up is summarized. 
Some of recent results at Nf = 2 are shown in Section |5[ The last section is devoted to a conclusion 
and outlook. 



2. Simulations at fixed topology 

2.1 Topology fixing term 

Dynamical simulations of the overlap fermions are quite nontrivial, because the sign function 



in Eq. (1.1) has a discontinuity at X = 0, where X is an eigenvalue of Hw- Vanishing X may 
occur during molecular-dynamics steps of the HMC update, at which the topological charge of 
the system changes its value. To keep the acceptance rate of HMC, one needs to take care of this 
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discontinuity so as to conserve the Hamiltonian precisely. One option is the reflection/refraction 
prescription [§]. This method first determines the time step at which X vanishes, and there the 
change of the pseudofermion action is computed. Like the light traveling across a surface of water, 
if the momentum of the mode is larger than the change of the pseudofermion action, it is refracted, 
while otherwise reflected. This method requires additional inversions of D(m) at A = 0, and thus 
the numerical cost quickly increases when the low-mode density becomes large. 

We instead employ a topology fixing term, which is implemented with an extra Wilson fermion 
and a twisted mass ghost as [^, [7|, ||] 

det ( ^2^2 ) = / ®Z*9Zov(-Se), (2-1) 



Se=X 



X- (2.2) 



This term is irrelevant in the continuum limit, and considered as a part of the gauge action. The 
numerator of Eq. ( |2.1| ) suppresses near-zero modes of Hw, while contribution from large frequency 
modes are compensated by the denominator. Since A = is prohibited, the topological charge 
is fixed during the molecular-dynamics update, and hence the reflection/refraction is no longer 
needed. We set the twisted ghost mass jj. = 0.2 throughout this work. As was shown in Refs. [f| 



10, 11 1, this term successfully suppresses the near-zero modes of Hw- 



As an alternative approach, the tunneling HMC has been proposed recently in Ref. [12]. This 
method also employs the extra Wilson fermion term, and thus near-zero modes of Hw are sup- 
pressed, but it projects out a few lowest-lying modes during the molecular-dynamics steps so as to 
enable topology changes. The tunneling HMC does not avoid A = 0, and the reflection/refraction 
prescription is necessary. Practical feasibility test is to be performed. 

2.2 Simulations at fixed topology 

Simulations at fixed topology are especially useful in the £ -regime where \-/m n y>> L (see 
Sec. |5T| ). In the £-regime, the topological charge dependence of physical observables is manifest 
and the fixed topology simulation is useful to determine the low-energy constants appearing in the 
effective chiral Lagrangian. On the other hand, in the ordinary regime (/^-regime), precision calcu- 
lations are possible only when the following two conditions are satisfied; (1) A physical observable 
in the 8 vacuum is related to those in the fixed-<2 vacua as a systematic expansion in terms of V -1 , 
where V is the spacetime volume. (2) The topological susceptibility, % t = (Q 2 )/V, is calculable and 
reproduces the known behavior from ChPT, which ensures that the local fluctuation of topological 
charge is active so as to produce relevant physics in a finite volume. Here we explain these points 



along Refs. [JlJ, |HJ [15p. 

What we want to compute is an expectation value in the 6 vacuum, which is related to the 
"vacua" with fixed values of Q through the Fourier transformation, 

Z(e)=X>-'' e %, Zfl = £g*»flZ(fl). (2.3) 

While the cluster property does hold for the former, the global topology becomes irrelevant in the 
infinite volume limit and the local fluctuation of topological charge is responsible to the physics. 
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In practice, however, the volume is inevitably finite. For a large enough volume, j;,V > 1, and for 
Q ^ XtV ■< the saddle point analysis is applicable and it leads to the conclusion that the distribution 
of Q is Gaussian and physical observables are represented as (0)q = (0}q + (finite V corrections). 
As an explicit example, a CP even correlator is represented as 



G Q = G(0) + G^(0)-^- 
XtV 



Q 2 c 4 



1 



XtV 2 X fV 



+ ^ (4) (0)^ + O(V- 3 ), (2.4) 



where is the n-th derivative of G(d) with respect to 6, and c 4 = -((Q A ) - 3{Q 2 ) 2 ) /V. To 

quantify the 0(1 /V) effect in Eq. ( |2.4j ), G^ must be known. This is determined with a help of the 
chiral perturbation theory (ChPT). For example, the 0(V~ l ) fixed topology effect to the PS meson 
mass is calculated in Ref. JT3|] at the tree-level of ChPT. It is also possible to determine G^- 1 ' in 



numerical simulations by comparing the results in different Q sectors. 

In the above formula, the topological susceptibility %t plays a key role. For a self-contained 
calculation, % t must be computed on the fixed-<2 configurations. % t is represented as a correlation of 



the local topological charge (o(x), which is also subject to Eq. ( |2.4[ ). From the clustering property 
in the 6 vacuum, 

UmJco(x)co(0)) Q = I - X t - + ^(V- 3 ). (2.5) 

Using the axial Ward-Takahashi identity, 

lim (mP°(x)mP°(0)) Q = lim (co(x)co(0)) Q , (2.6) 

\x\— >°° 

where P°(x) is the flavor singlet pseudoscalar density. Thus the topological susceptibility can be 
extracted from the long range behavior of the correlation function carrying the quantum number 
of rj' meson. If the determined % t in a simulation exhibits a reasonable value, it implies that in 
that system the local fluctuation of topological charge is active enough, and the system size is 



sufficiently large. Extraction of % t in our Nj = 2 simulation will be described in Sec. 5.2. Such 
a method is useful not only in the fixed topology simulations but also in standard HMC updates, 
because the changes of topological charge with a continuous variation of link variables become 
increasingly rare as approaching the continuum limit. 

These arguments indicate that the fixed topology simulations can provide a framework to de- 
termine the physical observables in the 8 vacuum in a self-contained manner with the finite size 
effects under control. 



3. Algorithm 

3.1 Overlap operator 

The overlap operator with a quark mass m is written as 

(in \ f m\ 
m ° + 2J + \ m °~ l) ^ si § n ^( _m °)]- ( 3 -^ 

mo is set to 1.6 throughout this work. The sign function means that an eigenmode (A, v^) of Hw is 
transformed to (sign (A), vi). Since the calculation of all the eigenmodes is impractical, some kind 
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of approximation of the sign function is required. We employ the Zolotarev rational approximation 

1 



[16, 17], 



X, 



(h W + C 2 n) }_ ■ 



,= lh w +C 2 l-l 



(3.2) 



where hyy = Hw /Knin with Xmin the eigenvalue having the smallest absolute value, do, c/, bi are 
easily calculable parameters. Although (hyp +c%i-\ ) _1 must be calculated N times, these terms can 



be obtained simultaneously by the multi-shift CG method [jig, 19]. Thus, the numerical cost mildly 
depends on N. This formula is valid in the region |A| € [A mi -„, X max ]. Since the error of the formula 
scales as exp(— X m i n N), the smaller A m( „ requires the larger N to keep the precision unchanged. 
Instead, if one calculates low-lying eigenvalues of Hw of \X\ < \h rs , one can determine the sign 
function of these modes explicitly and project them out from Hw- Then Xthrs replaces X m i„ in the 
above formula, leading to 



Nev 

sign(# w ) = Y, si g n (Ay) v i ® v j + sign(/fy)P ff , 



(3.3) 



where P 'h = 1 — Y!jZ\ v j ® v }> an d ^ev the number of modes with \Xj\ < Xthrs- The approximation 
formula is applied to the second term of this equation. The numerical cost depends on the density of 
low modes, since the determination of the low modes requires non-negligible computational time, 
and also projecting out the low modes from the sign function requires additional operations. Also in 
this sense, employing the topology fixing term improves the simulations of overlap fermions. In this 
work, we adopt Xthrs = 0.045 and Af = 10, which lead to an accuracy of |sign 2 //w — 1| — 10~( 7 ~ 8 ). 

3.2 Solver algorithms 

Since the overlap operator must be inverted at each step of the molecular-dynamics evolution, 
the improvement of the solver algorithm may significantly reduce the simulation cost. We have 
tested two algorithms; the nested CG (or 4DCG) method [|(| and the 5-dimensional CG (5DCG) 
||. The nested CG method is a straightforward implementation of the overlap solver. 



method [21 



It contains two nested CG iterations: an outer loop for operating D(m) and an inner loop for 
the calculation of {H w +C2/) _1 . The numerical cost can be reduced by applying the relaxation 
technique, which relaxes the convergence criterion of the inner loop as the outer loop iteration 
proceeds. Instead of the CG method, one can apply other Krylov subspace algorithms such as 
GMRES and SUMR. Since in HMC only the inversion of D f D appears and the CG method almost 
achieves the best performance, we compare the CG method with the 5D algorithm in the following. 

The 5-dimensional CG solver is based on the Schur decomposition [ ET] , [22| ]. Let us consider a 
5-dimensional block matrix (the N = 2 case is displayed as an example), 



Ms 











\ 




















Hw 















—Hw 









VP2 







R75+P0HJ 






(3.4) 
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Figure 1: Comparisons of the overlap solvers. The left panel shows the residual of the CG solver without 
the low-mode projection. The right panel shows the required time for the convergence in the case with the 
projection of 8 low modes of Hw- 



Since M5 can be decomposed as 

Ms=\ "_, , I I ' ' " II"'"" )■ (3.5) 




where S = D — CA B is called the Schur complement. One can solve a 4D linear equation S1//4 = 
Xa by solving a 5D equation 

-(:>(:)• 



Setting the parameters R, po, pi and (i = 1, . . . ,A0 in Eq. ( p.4| ) appropriately, the 5D solver can 
be used to invert the overlap operator approximated by Eq. (3.2). The 5D solver is accelerated 
by applying the even-odd preconditioning. Then one needs to solve a reduced linear equation, 
(1 -M~ e l M eo M m }M oe )\if e = x' e , where even/odd blocks of M 5 is denoted by M ee , M eo , etc. The 
inversions M~} and M^ 1 are easily calculated by forward/backward substitutions. The projection 
of low-modes is not straightforward for the even-odd preconditioned 5D solver, since the operation 
of M^ 1 (or MqJ) becomes much involved. Nevertheless it can be implemented cheaply, because 
the subspace of the matrix is spanned by x e , Jsx e , vj e , Y$Vj e (J = 1, .. . ,N ev ) [23]. 



Figure [T] compares the 4D and 5D solvers. The left panel shows the behavior of the residual 
against the numbers of the Wilson-Dirac operator multiplications, in the case of no low-mode 



projection [|28|]. The relaxation accelerates the nested CG by a factor of two, while the 5D solver 
exhibits much faster convergence for practical values of N. The required time to solve the linear 
equation is compared in the right panel for the case that the number of projected low modes is 8. 
This shows that the 5D solver is 3-4 times faster than the 4D solver in the whole region of quark 
mass used in this work. Therefore we mainly use the 5D solver in HMC. For the computation of the 
quark propagator, we adopt the nested 4D solver because of an advantage to obtain the propagators 
with several valence quark masses simultaneously with the multi-shift CG algorithm. 

Here let us compare the cost of simulations with the domain-wall fermion. Let us consider the 
cost to solve a linear equation Dx = b. At m ~ m s /2 ,aT x ~ 1.7 GeV, and on 16 3 x 32 lattices, the 
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domain-wall solver requires 0(800) iterations for a precision \r\/\b\ < 10~ 10 according to Ref. [|24|], 
while the overlap solver requires 0(1200) iterations. To convert to the numbers of the Wilson- 
Dirac operator multiplications, 2N S and 2N + 1 are multiplied for the domain-wall and the overlap 
fermions, respectively. Then the difference amounts to a factor of about 2.5. In HMC, the overlap 
solver must be called twice, while the domain-wall solver is applied to D^D. This leads to another 
factor of 2 difference leading to the total difference of 0(5). 

3.3 Hybrid Monte Carlo algorithm 

Implementation of the HMC algorithm for the overlap operator with the approximation for- 



mula (3.2) is straightforward except for the discontinuity at A = 0. Since we employ the topology 



fixing term, A = does not appear and no reflection/refraction prescription is required. In order 



to improve the performance of HMC, we adopt the mass preconditioning [25] together with the 
multi-time step procedure [|26|]. Introducing a preconditioning term with a heavier quark mass m' 
than the dynamical quark mass m, the fermion action becomes S F = Sp F i + Sp F 2, 

S PF1 = ^[D{m')^D(m')]-%, S PF2 = $ {D(m')[D(m)^ D{m)}- l D(m^} <fc. (3.7) 

The forces from the preconditioner (PF1), the preconditioned dynamical quark (PF2), the gauge 
field (G), and the extra Wilson fermion/ghost field (E) have a hierarchical structure, 

Fq ~ Fp S> Fpp\ 3> Fpp2- (3-8) 

Thus we set 

At (PF2) > Ar(p n ) » At (g) = Ax {E) . (3.9) 

AT(£) is set to be the same as the gauge part also to ensure the disappearance of near-zero modes 
of H w . The cost to calculate F E is negligible compared to the overlap fermions. In the case 
of Nf = 2 + 1, the time steps of the two-flavor and one-flavor parts are set equal. In this work, 
At(p F2 )/AT(p/t i) and Ax/ PF1 \ / Ax/ G ^ are set to 4-6. 



In the Nf = 2 simulation, we employ the noisy Metropolis test [|27J together with a less precise 



5D solver without the projection of low-modes of Hw [10, 28]. During the molecular dynamics, 
the calculation of the low-lying modes of Hw is skipped. We fix the value of A,/ irv , the lower-bound 
of the region where the approximation formula ( |3T2| ) is valid, though the modes with \X\ < Xth rs 
may appear. The error with this setting is corrected by the noisy Metropolis test performed at the 
end of each MD evolution in addition to the usual Metropolis test. This algorithm is twice faster 
than the case with the 4D solver with the projection of low-modes of Hw, which was used in an 
early stage of the simulation. 



The Nf = 2 + 1 simulation has been started recently []23[]. The one-flavor part is implemented 



with one of the chirality sectors [g^, [30p making use of the fact that H(m) commutes with 75, thus 



H 2 =P+H 2 P + +P H 2 P =Q+ + Q , detH 2 = det 0+ • det 0_ . (3.10) 
Except for the trivial contribution from the zero-modes, the determinant of one chirality sector cor- 

4- 1 

responds to the contribution of one flavor. Thus, the pseudofermion action S F = <j>oQa <Po> where 
a can be either + or — , represents the one-flavor of dynamical fermion. The same acceleration 
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« « I 1 

I 
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in 



J 
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Nf=0, Iwasaki (CP-PACS, 2004) 
Nf=0. Iwasaki + extra Wilson 
Nf=2, Iwasaki + ex Wilson + overlap 
Nf=2, Iwasaki + clover (CP-PACS, 2002) 



Figure 2: The left panel shows the lattice scale a(ro) set by ro = 0.49fm. The right panel shows the j3- 
dependence of a(ro) for Nf = 2 and 0, together with the result of the clover fermion. 



techniques as Nf = 2 are applicable to the one-flavor part. For the Nf = 2 + 1 simulation, we 
adopt the 5D solver with the projection of low-modes of Hw, which no longer requires the noisy 



Metropolis test [23]. 



4. Simulation 



Numerical simulations are performed mainly on IBM Blue Gene/L system at KEK. At present 
the sustained performance on one rack of Blue Gene/L (1024 nodes, 5.7 TFlops of peak perfor- 



mance) is about 30% for the Wilson operator [31 1, and 10-15% for the overlap HMC. 

The Nf = 2 simulations are performed on 16 3 x 32 lattices at j8 = 2.3. We use 6 quark masses, 



0.015, 0.025, 0.035, 0.050, 0.070, and 0.100, roughly corresponding to mf ys /6 
trajectories of a length 0.5 are generated at each sea quark mass in the Q = sector. For m uc i = 0.050 
(~ mf /2), configurations in Q = — 2 and —4 sectors are also generated. In addition, we perform 
a simulation in the £-regime at j3 = 2.35 with a quite small quark mass, m uc i ~ 3 MeV. To generate 
one trajectory, 1 1-26 minutes are needed on one rack of Blue Gene/L. The acceptance rate is kept 



mf ys . 10,000 



to 80-90% at each quark mass. The locality of the overlap operator was examined in Ref. [|32[]. 

The Nf = 2 + 1 simulations have been started recently on 16 3 x 48 lattices with almost the same 
parameters as Nf = 2, while each trajectory has a length 1.0 [23]. We use 5 values of m u d covering 
the same quark mass region as the Nf = 2 case for each of 2 strange quark masses, m s = 0.080 and 
0.100, around the physical strange quark mass m p s hys . Present performance is around 2 hours for 
one trajectory on one rack of Blue Gene/L. 

The lattice scale is set by the hadronic radius ro which is defined through 



,V(r) 



1.65, 



(4.1) 



by setting the physical value ro = 0.49 fm. The static quark potential V(r) is calculated with 
the standard procedure. Figure || shows the result of the lattice spacing. The left panel displays 
the result for Nf = 2 extrapolated to the chiral limit, as well as a preliminary result for Nf = 
2 + 1. A linear extrapolation of the Nf = 2 data gives a(m uc i = 0) = 0.1184(12)^(11)^ fm. 
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The right panel shows the /? -dependence of a in the chiral limit for Nf = and 2 together with 
the clover fermion case. The shift of /3 (a) with respect to the number of flavors is milder for the 
overlap fermion than for the Wilson-type fermions. This behavior is consistent with the perturbative 
calculation. 



5. Results 

On the Nf = 2 lattices, we finished the generation of gauge configurations and are currently 
calculating various physical observables. The following calculations are in progress. 
• £ -regime 



• Topological Susceptibility [37, p8|] 

• Pion mass and decay constant [ p9[ | 

• Pion form factor p0| ] 
B meson bag parameter p|] 
Tl + -K Q mass difference [P 



Pion scattering length [ 43 ] 



In the following, we briefly describe the first three subjects, which are the first testing ground of 
the viability of our simulations. 

5.1 £ -regime 

The chiral condensate is related to the spectral density of the Dirac operator through the Banks- 



Casher relation Q44[], 

L = -(^)= lim lim^p, (5.1) 

where p(A) = Jjt(5(A — is the spectral density of the Dirac operator. The accumulation of 
low modes generate the spontaneous chiral symmetry breaking. Two limits, m — ► oo and V — > °°, 
in the above equation are not commutable. However, it is also convenient to consider the opposite 
order of the limit. The £ -regime is defined through the condition 

l/A QCD <L<l/m^ (5.2) 

which implies m -C l/LV. In the £-regime, the low energy effective theory is applied with the 
same parameters as the infinite V case. Because of the finiteness of the volume, the topological 
charge dependence of observables becomes manifest, and hence it is convenient to determine the 
parameters of the low energy effective Lagrangian in the fixed topology simulations. Another 
advantage of the £-regime simulations is that the chiral random matrix theory (RMT) is expected 
to describe the behavior of the low-lying modes. 

The results in the £-regime for Nf = 2 are displayed in Figure [3] [|34|, 35]. We find good 



agreement with the RMT prediction for the level distribution. The lowest level distribution gives a 
value of the chiral condensate which is consistent with other determinations. Another agreement 
with the RMT prediction is found for the so-called topology-flavor duality, as shown in the right 
panel of Figure [| With the nonperturbative renormalization, we obtaine 

I M (2GeV) = (251±7(stat)± ll(syst) MeV) 3 . (5.3) 

The main source of the systematic error is <9(£ 2 ), which can be corrected using meson correlators 

ph. 
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10 20 30 40 50 60 Q 

VLV 



Figure 3: The result in the e-regime atiVV = 2. The left panel shows the spectral density of the overlap-Dirac 
operator. The right panel compares the level spacings in the cases of different Q and Nf. 



5.2 Topological susceptibility 

Now let us go back to the /^-regime. We first need to show that the topological susceptibility 



is successfully determined in the fixed-<2 "vacua" from the correlation function (Sec. 2.2). The 
topological susceptibility Xt is extracted from the correlation function 



1 



C 11 >(t) = JJ J £(m (j P (x)m lj P\0)) Q , 

X 

where P°(x) is the flavor singlet pseudoscalar density. For the overlap fermions, 



(5.4) 



1 



Nf 



p0 « = w I 



aD{m q = 0) 
2m 



At large t, 



<V(0 



1 

V 



V 



Xt 



2XtV 



+ 



-3 



) + ■ 



Y f (x). 



(e- m n") 



(5.5) 



(5.6) 



hence the value of % t is extracted from the plateau of CV(f). 

The measurement of % t is performed on Nf = 2 and Q = lattices at every 20 trajectories, thus 
with statistics of 500. Prior to measuring the correlator, we calculate 50 pairs of low-lying eigen- 
modes of the overlap operator D(0) by the implicitly restarted Lanczos algorithm. Since these 
low-modes can be explicitly inverted, the solver algorithm is applied to the operator projected out 
these modes, whose condition number is reduced from the original one. This low-mode precon- 
ditioning accelerates the calculation of the quark propagators by a factor of 8. These low-modes 
are also used in the low-mode averaging, which is an average of correlators over all the spacetime 
source points for the low-mode contributions p5| , fff]]. For the disconnected part of the correlator 



(p.4|), the quark propagator is approximated by 50 pairs of the low-eigenmodes, by observing these 
modes dominate the correlator. 

The result is shown in Figure^. The left panel displays the correlator — C^i(t) at m q = 0.025. 



The data are fitted to a function A +B(e +e 



Assuming |c4| <C 2% t V , we obtain a Xt 



3.40(27) x 10 at m q = 0.025. The right panel shows the extracted % t as a function of sea quark 
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Figure 4: The result of the topological susceptibility for the Nf = 2 simulation. The left panel shows the 
extraction of the topological susceptibility %, from the 77' correlator at m q = 0.025. The right panel displays 
% t rt versus the sea quark mass m q r§. 



mass in units of rn. For the smallest three quark masses, the data are well fitted to a linear function 
whose intercept is consistent with zero. This behavior is consistent with the prediction of the chiral 
perturbation theory [|47|], 

Xt = ^ + 0{m 2 q ). (5.7) 

From the slope of the fit result, a value of the chiral condensate is obtained as r^L = 0.182(6). In 
order to convert it to a physical value, we use the renormalization factor Z ; ^ s (2GeV) = 0.742(12), 
which is obtained nonperturbatively through the RI/MOM scheme on the lattice |39|]. This 
leads to a value 

£ M (2GeV) = (254±5(stat)±10(syst)MeV) 3 , (5.8) 



which is consistent with the value in Eq. (5.3) independently obtained in the £-regime. The statisti- 
cal error includes those of aT x and Z^f s , and the systematic error is of the higher order effects such 
as the C4 term. 

We also measure % t on the Nf = 2 lattices with Q = —2 and —4 at m q = 0.050 with statistics 
of 250 configurations. The extracted values of %t w& consistent with the value at Q = 0. 

These results indicate that the local fluctuation of the topological charge is active enough to 
produce relevant chiral dynamics on the fixed-Q vacua. The successful extraction of Xt enables us 



to quantify the finite size effects for other observables represented in Eq. (2.4). 



5.3 Pion mass and decay constant 

Among hadronic observables, the pion mass and decay constant are the primary quantities to 
be tested with the chiral perturbation theory. Here we briefly present our results of the Nf = 2 



simulation for m % and f n [39]. 
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The meson correlators are computed at every 20 trajectories. The low-mode preconditioning 
and the low-mode averaging described in the previous subsection are applied to the meson correla- 
tors. The pion mass and decay constant are extracted from the pseudoscalar meson correlators with 
point and smeared sources by a simultaneous fit. f K is obtained through the axial Ward-Takahashi 
identity, 

f„ = 2m q (0\P(0)\7Z)/ml, (5.9) 

without further renormalization. Note that we are using the f n = 130 MeV normalization. The 
quark mass is renormalized with the renormalization factor Z^ s (2GeV) = 0.742(12), which is 



obtained nonperturbatively through the RI/MOM scheme on the lattice [J48|, p9p . 

The pion mass and decay constant receive two kinds of finite size effect (FSE): the standard 
FSE and that from fixed Q. The former effect is evaluated by the formulae determined by Colan- 
gelo et al. [fi9|], which is developed from the Liischer's formula for the relation between the pion 



scattering amplitude and the mass-shift in a finite box J50|]. We use the NNLO results of Ref. [ |49j ] 
and the values of low energy constants estimated in Ref. [51]. For our lightest quark mass, they 
amount to 4.5% and 6.0% for m\jm q and f n , respectively. The FSE from the fixed topological 



charge is subject to Eq. (2.4). The 0(V ) corrections to m n and f n are determined at NLO of 



ChPT The correction to m n starts at the tree-level. In the small quark mass region, it shifts the 
value of m| / m q with the same size as the standard FSE while in the opposite direction. Thus the 
two effects almost cancel each other. The fixed-<2 correction to f K starts at NLO of ChPT, and gives 
small effect on the result of f n . 

Figure || shows the results of nr\/m q and f K after correcting FSEs. To compare with the ChPT 
expansion, we chose £, = (m n / 4n f n ) 2 as an expansion parameter, where f n is the measured value 
at each quark mass. We apply three types of ChPT fit: NLO, NNLO, and a simplified NNLO 
called NNLO' hereafter. The NNLO' formulae are obtained from those of NNLO by using an 
approximation <^ 2 ln<^ — 2.5^ 2 , which is numerically reasonable in our target range, < <^ < 0.1. 
Explicit functional forms are given in Ref. |[39|]. The results of fits are displayed in Figure We 
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Figure 6: Comparisons of physical quantities obtained from the chiral fit ansatze: / (upper left), E 1 / 3 
(upper right), /? hys (lower left) and l^ hyli (lower right). In each panel, the red circles and blue squares are 
corresponding to the fits with 6 and 5 lightest data points, respectively. 



note that the fit to the NNLO formulae is performed by a simultaneous fit for m 2 n jm q and f n , and 
thus the value of % 2 cannot be compared directly to other fits. 

These fits determine the following quantities: f, L = Bq ■ f 2 /2, Z| hys , and Zj hys . Figure || 
compares the results of fits with the phenomenologically determined values (/ [52] and /| hys Z^ hys 



[pl[]). The result of £ is compared with our calculation in the £-regime in Sec. |5J] [p4[]. Except 
for the case of ys , the results of NLO fit are inconsistent with NNLO and the phenomenological 
estimates. This implies the failure of the NLO formulae to describe our data. The results of NNLO 
and NNLO' fits are consistent with each other and with the phenomenological estimates. 

To quantify the FSEs more accurately, we need simulations with a larger lattice size, as well 
as a comparison with the result at Q ^ 0. We emphasize that these results indicate that the overlap 
simulations with fixed topology provide a framework for precision calculations of the spectrum and 
the matrix elements in the chiral regime. 



6. Conclusion 

We are performing large-scale simulations with Nf = 2 and Nf = 2 + 1 dynamical overlap 
fermions. The range of sea quark mass covers m p s hys /6 - m% hys . The simulations are performed in 
fixed topological charge sectors. The results of the topological susceptibility and the pion mass and 
decay constant indicate that these simulations can provide a ground for precision computations of 
matrix elements with controlled chiral extrapolation. 

We completed a generation of Nf = 2 gauge configurations on 16 3 x 32 lattices with a ~ 0.12 
fm. Numbers of measurements are in progress and planned. The Nf = 2+1 simulations on 1 6 3 x 48 
lattices are also in progress. These configurations will be supplied to ILDG soon after the first 
publication of the result of the spectrum. 

For further investigation of finite size effects and for more extended objects than mesons, 
simulations at larger lattice sizes are desired. The target size of the spatial extent is 24, which 
requires further improvements of numerical algorithms. 
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